\defmodule{DigitalNet}

This class provides the basic structures for storing and
manipulating \emph{linear digital nets in base $b$}, for an arbitrary
base $b\ge 2$.  We recall that a net contains $n = b^k$ points in
$s$ dimensions, where the $i$th point $\mathbf{u}_i$,
for $i=0,\dots,b^k-1$, is defined as follows:
\begin{latexonly}
\begin{eqnarray*}
  i &=& \sum_{\ell=0}^{k-1} a_{i,\ell} b^\ell, \\ %  \eqlabel{eq:digital-i} \\
 \pmatrix{u_{i,j,1}\cr u_{i,j,2}\cr \vdots \cr }
    &=& \bC_j \pmatrix{a_{i,0}\cr a_{i,1}\cr \vdots \cr a_{i,k-1}\cr},
                                           \\ %  \eqlabel{eq:digital-Cj} \\
 u_{i,j} &=& \sum_{\ell=1}^\infty u_{i,j,\ell} b^{-\ell},
                                 \\[7pt] %  \eqlabel{eq:digital-uij} \\
  \bu_i &=& (u_{i,0},\dots,u_{i,s-1}).     \\ %  \eqlabel{eq:digital-ui}
\end{eqnarray*}
\end{latexonly}
\begin{htmlonly}
\begin{eqnarray}
  i &\quad=\quad& \sum_{r=0}^{k-1} a_{i,r} b^r, \\
  (u_{i,j,1}\  u_{i,j,2}\  \ldots )^T
    &\quad=\quad& \mathbf{C}_j\ (a_{i,0}\  a_{i,1}\  \ldots \ a_{i,k-1})^T, \\
 u_{i,j} &\quad=\quad& \sum_{r=1}^\infty u_{i,j,r} b^{-r},  \\
  \mathbf{u}_i &\quad=\quad& (u_{i,0},\ \ldots,\ u_{i,s-1}).
\end{eqnarray}
\end{htmlonly}
In our implementation, the matrices $\mathbf{C}_j$ are $r\times k$,
so the expansion of $u_{i,j}$ is truncated to its first $r$ terms.
The points are stored implicitly by storing the generator matrices
$\mathbf{C}_j$ in a large two-dimensional array of integers,
 with $srk$ elements.
\begin{detailed}
For general $b$, the element $(l,c)$ of $\mathbf{C}_j$ (counting elements
 from 0)  is stored at position $[jk+c][l]$ in this array.
\end{detailed}

The points $\mathbf{u}_i$ are enumerated using the Gray code
technique\latex{ as proposed in \cite{rANT79a,rTEZ95a}
(see also \cite{fGLA04a,vHON03a})}.
With this technique, the $b$-ary representation of $i$,
$\mathbf{a}_i = (a_{i,0}, \dots, a_{i,k-1})$, is replaced \latex{in Equation
(\ref{eq:digital-Cj})} by a Gray code representation of $i$,
$\mathbf{g}_i = (g_{i,0}, \dots, g_{i,k-1})$. The Gray code $\mathbf{g}_i$
 used here is defined by $g_{i,k-1} = a_{i,k-1}$ and
\latex{
$g_{i,\ell} = (a_{i,\ell} - a_{i,\ell+1}) \mod b$ for $\ell = 0,\dots,k-2$.}
\html{$g_{i,\nu} = (a_{i,\nu} - a_{i,\nu+1}) \mod b$ for $\nu = 0,\dots,k-2$.}
It has the property that
$\mathbf{g}_i = (g_{i,0}, \dots, g_{i,k-1})$ and
$\mathbf{g}_{i+1} = (g_{i+1,0}, \dots, g_{i+1,k-1})$
differ only in the position of the smallest index
\latex{$\ell$}\html{$\nu$} such that
\latex{$a_{i,\ell} < b - 1$, and we have $g_{i+1,\ell} = (g_{i,\ell}+1) \mod b$}\html{$a_{i,\nu} < b - 1$, and we have $g_{i+1,\nu} = (g_{i,\nu}+1) \mod b$}
in that position.

This Gray code representation permits a more efficient enumeration
of the points by the iterators.
It changes the order in which the points $\mathbf{u}_i$ are enumerated,
but the first $b^m$ points remain the same for every integer $m$.
The $i$th point of the sequence with the Gray enumeration
is the $i'$th point of the original enumeration, where $i'$ is the
integer whose $b$-ary representation $\mathbf{a}_{i'}$ is given by the Gray
code $\mathbf{g}_i$.
To enumerate all the points successively, we never need to compute
the Gray codes explicitly.
It suffices to know the position \latex{$\ell$}\html{$\nu$} of the Gray code digit
that changes at each step, and this can be found quickly from the
$b$-ary representation $\mathbf{a}_i$.
The digits of each coordinate $j$ of the current point can be updated by
adding column \latex{$\ell$}\html{$\nu$} of the generator matrix $\mathbf{C}_j$ to the old digits,
modulo $b$.

One should avoid using the method \method{getCoordinate}{}{(i, j)}
for arbitrary values of \texttt{i} and \texttt{j}, because this is much
slower than using an iterator to access successive coordinates.

Digital nets can be randomized in various ways
\cite{mMAT99a,rFAU02a,vLEC02a,vOWE03a}.
Several types of randomizations specialized for nets are implemented
directly in this class.

A simple but important randomization is the \emph{random digital shift}
in base $b$, defined as follows: replace each digit $u_{i,j,\latex{\ell}\html{\nu}}$ in \html{the third equation above}
\latex{(\ref{eq:digital-uij})} by $(u_{i,j,\latex{\ell}\html{\nu}} + d_{j,\latex{\ell}\html{\nu}}) \mod b$,
where the $d_{j,\latex{\ell}\html{\nu}}$'s are i.i.d.\ uniform over $\{0,\dots,b-1\}$.
This is equivalent to applying a single random shift to all the points
in a formal series representation of their coordinates \cite{vLEC02a,vLEM03a}.
In practice, the digital shift is truncated to $w$ digits,
for some integer $w\ge r$.
Applying a digital shift does not change the equidistribution
and $(t,m,s)$-net properties of a point set \cite{vHON03a,vLEC99a,vLEM03a}.
Moreover, with the random shift, each point is uniformly distributed over
the unit hypercube (but the points are not independent, of course).

A second class of randomizations specialized for digital nets
are the \emph{linear matrix scrambles}
\cite{mMAT99a,rFAU02a,vHON03a,vOWE03a}, which multiply the matrices
 $\mathbf{C}_j$
by a random invertible matrix $\mathbf{M}_j$, modulo $b$.
There are several variants, depending on how $\mathbf{M}_j$ is generated, and on
whether $\mathbf{C}_j$ is multiplied on the left or on the right.
In our implementation, the linear matrix scrambles are incorporated directly
into the matrices $\mathbf{C}_j$\latex{ (as in \cite{vHON03a})}, so they
 do not slow down  the enumeration of points.
Methods are available for applying
linear matrix scrambles and for removing these randomizations.
These methods generate the appropriate random numbers and make
the corresponding changes to the $\mathbf{C}_j$'s.
A copy of the original $\mathbf{C}_j$'s is maintained, so the point
set can be returned to its original unscrambled state at any time.
When a new linear matrix scramble is applied, it is always applied to
the \emph{original} generator matrices.
% (so the new matrix scramble is combined with the previous ones).
The method \method{resetGeneratorMatrices}{} removes the current matrix
scramble by resetting the generator matrices to their original state.
On the other hand, the method \method{eraseOriginalGeneratorMatrices}{}
replaces the original generator matrices by the current
ones, making the changes permanent.
This is useful if one wishes to apply two or more linear matrix
scrambles on top of each other.

Linear matrix scrambles are usually combined with a random digital shift;
this combination is called an \emph{affine matrix scramble} \cite{vOWE03a}.
These two randomizations are applied via separate methods.
The linear matrix scrambles are incorporated into the matrices $\mathbf{C}_j$
whereas the digital random shift is stored and applied separately,
independently of the other scramblings.

Applying a digital shift or a linear matrix scramble to a digital net
invalidates all iterators for that randomized point,
because each iterator uses a
\emph{cached} copy of the current point, which is updated only when
the current point index of that iterator changes, and the update also
depends on the cached copy of the previous point.
After applying any kind of scrambling, the iterators must be
reinitialized to the \emph{initial point} by invoking
\externalmethod{}{PointSetIterator}{resetCurPointIndex}{}
or reinstantiated by the \method{iterator}{} method
(this is not done automatically).


\bigskip\hrule\bigskip
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{code}
\begin{hide}
/*
 * Class:        DigitalNet
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.hups;\begin{hide}

import umontreal.iro.lecuyer.util.PrintfFormat;
import umontreal.iro.lecuyer.rng.*;
\end{hide}

public class DigitalNet extends PointSet \begin{hide} {

   // Variables to be initialized by the constructor subclasses.
   protected int b = 0;         // Base.
   protected int numCols = 0;   // The number of columns in each C_j. (= k)
   protected int numRows = 0;   // The number of rows in each C_j. (= r)
   protected int outDigits = 0; // Number of output digits (= w)
   private int[][] originalMat; // Original gen. matrices without randomizat.
   protected int[][] genMat;    // The current generator matrices.
                                // genMat[j*numCols+c][l] contains column c
                                // and row l of C_j.
   protected int[][] digitalShift; // The digital shift, initially zero (null).
                                // Entry [j][l] is for dimension j, digit l,
                                // for 0 <= l < outDigits.
   protected double normFactor; // To convert output to (0,1); 1/b^outDigits.
   protected double[] factor;   // Lookup table in ascending order: factor[i]
                                // = 1/b^{i+1} for 0 <= i < outDigits.

   // primes gives the first index in array FaureFactor
   // for the prime p. If primes[i] = p, then
   // FaureFactor[p][j] contains the Faure ordered factors of base p.
   private int[] primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,
      43, 47, 53, 59, 61, 67};

   // Factors on the diagonal corresponding to base b = prime[i] ordered by
   //  increasing Bounds.
   private int[][] FaureFactor = {{1}, {1, 2}, {2, 3, 1, 4},
      { 2, 3, 4, 5, 1, 6}, {3, 4, 7, 8, 2, 5, 6, 9, 1, 10},
      { 5, 8, 3, 4, 9, 10, 2, 6, 7, 11, 1, 12},
      { 5, 7, 10, 12, 3, 6, 11, 14, 4, 13, 2, 8, 9, 15, 1, 16},
      { 7, 8, 11, 12, 4, 5, 14, 15, 3, 6, 13, 16, 2, 9, 10, 17, 1, 18},
      { 5, 9, 14, 18, 7, 10, 13, 16, 4, 6, 17, 19, 3, 8, 15, 20, 2, 11, 12,
        21, 1, 22},
      { 8, 11, 18, 21, 12, 17, 9, 13, 16, 20, 5, 6, 23, 24, 4, 7, 22, 25, 3,
        10, 19, 26, 2, 14, 15, 27, 1, 28},
      { 12, 13, 18, 19, 11, 14, 17, 20, 7, 9, 22, 24, 4, 8, 23, 27, 5, 6, 25,
        26, 3, 10, 21, 28, 2, 15, 16, 29, 1, 30},
      { 8, 14, 23, 29, 10, 11, 26, 27, 13, 17, 20, 24, 7, 16, 21, 30, 5, 15,
        22, 32, 6, 31, 4, 9, 28, 33, 3, 12, 25, 34, 2, 18, 19, 35, 1, 36},
      { 16, 18, 23, 25, 11, 15, 26, 30, 12, 17, 24, 29, 9, 32, 13, 19, 22,
        28, 6, 7, 34, 35, 5, 8, 33, 36, 4, 10, 31, 37, 3, 14, 27, 38, 2, 20,
        21, 39, 1, 40},
      { 12, 18, 25, 31, 9, 19, 24, 34, 8, 16, 27, 35, 10, 13, 30, 33, 15, 20,
        23, 28, 5, 17, 26, 38, 6, 7, 36, 37, 4, 11, 32, 39, 3, 14, 29, 40, 2,
        21, 22, 41, 1, 42},
      { 13, 18, 29, 34, 11, 17, 30, 36, 10, 14, 33, 37, 7, 20, 27, 40, 9, 21,
        26, 38, 15, 22, 25, 32, 6, 8, 39, 41, 5, 19, 28, 42, 4, 12, 35, 43,
        3, 16, 31, 44, 2, 23, 24, 45, 1, 46},
      { 14, 19, 34, 39, 23, 30, 12, 22, 31, 41, 8, 11, 20, 24, 29, 33, 42, 45,
        10, 16, 37, 43, 7, 15, 38, 46, 17, 25, 28, 36, 5, 21, 32, 48, 6, 9,
        44, 47, 4, 13, 40, 49, 3, 18, 35, 50, 2, 26, 27, 51, 1, 52},
      { 25, 26, 33, 34, 18, 23, 36, 41, 14, 21, 38, 45, 24, 27, 32, 35, 11,
        16, 43, 48, 9, 13, 46, 50, 8, 22, 37, 51, 7, 17, 42, 52, 19, 28, 31,
        40, 6, 10, 49, 53, 5, 12, 47, 54, 4, 15, 44, 55, 3, 20, 39, 56, 2,
        29, 30, 57, 1, 58},
      { 22, 25, 36, 39, 17, 18, 43, 44, 24, 28, 33, 37, 13, 14, 47, 48, 16,
        19, 42, 45, 9, 27, 34, 52, 8, 23, 38, 53, 11, 50, 7, 26, 35, 54, 21,
        29, 32, 40, 6, 10, 51, 55, 5, 12, 49, 56, 4, 15, 46, 57, 3, 20, 41,
        58, 2, 30, 31, 59, 1, 60},
      { 18, 26, 41, 49, 14, 24, 43, 53, 12, 28, 39, 55, 29, 30, 37, 38, 10,
        20, 47, 57, 16, 21, 46, 51, 8, 25, 42, 59, 13, 31, 36, 54, 9, 15,
        52, 58, 7, 19, 48, 60, 23, 32, 35, 44, 5, 27, 40, 62, 6, 11, 56,
        61, 4, 17, 50, 63, 3, 22, 45, 64, 2, 33, 34, 65, 1, 66}
};


   public double getCoordinate (int i, int j) {
      // convert i to Gray representation, put digits in gdigit[].
      int l, c, sum;
      int[] bdigit = new int[numCols];
      int[] gdigit = new int[numCols];
      int idigits = intToDigitsGray (b, i, numCols, bdigit, gdigit);
      double result = 0;
      if (digitalShift != null && dimShift < j)
         addRandomShift (dimShift, j, shiftStream);
      for (l = 0; l < outDigits; l++) {
         if (digitalShift == null)
            sum = 0;
         else
            sum = digitalShift[j][l];
         if (l < numRows)
            for (c = 0; c < idigits; c++)
               sum += genMat[j*numCols+c][l] * gdigit[c];
         result += (sum % b) * factor[l];
      }
      if (digitalShift != null)
         result += EpsilonHalf;
      return result;
   }

   public PointSetIterator iterator() {
      return new DigitalNetIterator();
   }\end{hide}
\end{code}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructor}
\begin{code}

   public DigitalNet ()\begin{hide} {
     }\end{hide}
\end{code}
\begin{tabb} Empty constructor.
\end{tabb}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}
\begin{code}

   public double getCoordinateNoGray (int i, int j)\begin{hide} {
      // convert i to b-ary representation, put digits in bdigit[].
      int l, c, sum;
      int[] bdigit = new int[numCols];
      int idigits = 0;
      for (c = 0; i > 0; c++) {
         idigits++;
         bdigit[c] = i % b;
         i = i / b;
      }
      if (digitalShift != null && dimShift < j)
         addRandomShift (dimShift, j, shiftStream);
      double result = 0;
      for (l = 0; l < outDigits; l++) {
         if (digitalShift == null)
            sum = 0;
         else
            sum = digitalShift[j][l];
         if (l < numRows)
            for (c = 0; c < idigits; c++)
               sum += genMat[j*numCols+c][l] * bdigit[c];
         result += (sum % b) * factor[l];
      }
      if (digitalShift != null)
         result += EpsilonHalf;
      return result;
   }\end{hide}
\end{code}
 \begin{tabb}
   Returns $u_{i,j}$, the coordinate $j$ of point $i$, the points
   being enumerated in the standard order (no Gray code).
 \end{tabb}
\begin{htmlonly}
   \param{i}{point index}
   \param{j}{coordinate index}
   \return{the value of $u_{i,j}$}
\end{htmlonly}
\begin{code}

   public PointSetIterator iteratorNoGray()\begin{hide} {
      return new DigitalNetIteratorNoGray();
   }\end{hide}
\end{code}
\begin{tabb}
  This iterator does not use the Gray code. Thus the points are enumerated
  in the order of their first coordinate before randomization.
\end{tabb}
\begin{code}

   public void addRandomShift (int d1, int d2, RandomStream stream)\begin{hide} {
      if (null == stream)
         throw new IllegalArgumentException (
              PrintfFormat.NEWLINE +
              "   Calling addRandomShift with null stream");
      if (0 == d2)
         d2 = Math.max (1, dim);
      if (digitalShift == null) {
         digitalShift = new int[d2][outDigits];
         capacityShift = d2;
      } else if (d2 > capacityShift) {
         int d3 = Math.max (4, capacityShift);
         while (d2 > d3)
            d3 *= 2;
         int[][] temp = new int[d3][outDigits];
         capacityShift = d3;
         for (int i = 0; i < d1; i++)
            for (int j = 0; j < outDigits; j++)
               temp[i][j] = digitalShift[i][j];
         digitalShift = temp;
      }
      for (int i = d1; i < d2; i++)
         for (int j = 0; j < outDigits; j++)
            digitalShift[i][j] = stream.nextInt (0, b - 1);
      dimShift = d2;
      shiftStream = stream;
     }\end{hide}
\end{code}
\begin{tabb}  Adds a random digital shift to all the points of the point set,
  using stream \texttt{stream} to generate the random numbers.
  For each coordinate $j$ from \texttt{d1} to \texttt{d2-1},
  the shift vector $(d_{j,0},\dots,d_{j,k-1})$
  is generated uniformly over $\{0,\dots,b-1\}^k$ and added modulo $b$ to
  the digits of all the points.
  After adding a digital shift, all iterators must be reconstructed or
  reset to zero.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used to generate uniforms}
\end{htmlonly}
\begin{code}

   public void addRandomShift (RandomStream stream)\begin{hide} {
      addRandomShift (0, dim, stream);
     }\end{hide}
\end{code}
\begin{tabb}  Same as \method{addRandomShift}{}\texttt{(0, dim, stream)},
  where \texttt{dim} is the dimension of the digital net.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used to generate uniforms}
\end{htmlonly}
\begin{hide}\begin{code}

   public void clearRandomShift() {
      super.clearRandomShift();
      digitalShift = null;
   }
\end{code}
\begin{tabb}
   Erases the current digital random shift, if any.
\end{tabb}\end{hide}
\begin{code}\begin{hide}

   public String toString() {
      StringBuffer sb = new StringBuffer (100);
      if (b > 0) {
         sb.append ("Base = ");   sb.append (b);
         sb.append (PrintfFormat.NEWLINE);
      }
      sb.append ("Num cols = ");   sb.append (numCols);
      sb.append (PrintfFormat.NEWLINE + "Num rows = ");
      sb.append (numRows);
      sb.append (PrintfFormat.NEWLINE + "outDigits = ");
      sb.append (outDigits);
      return sb.toString();
   }


   // Print matrices M for dimensions 0 to N-1.
   private void printMat (int N, int[][][] A, String name) {
      for (int i = 0; i < N; i++) {
         System.out.println ("-------------------------------------" +
            PrintfFormat.NEWLINE + name + "   dim = " + i);
         int l, c;   // row l, column c, dimension i for A[i][l][c].
         for (l = 0; l < numRows; l++) {
            for (c = 0; c < numCols; c++) {
               System.out.print (A[i][l][c] + "  ");
            }
            System.out.println ("");
         }
      }
      System.out.println ("");
   }


   // Print matrix M
   private void printMat0 (int[][] A, String name) {
         System.out.println ("-------------------------------------" +
                             PrintfFormat.NEWLINE + name);
         int l, c;   // row l, column c for A[l][c].
         for (l = 0; l < numCols; l++) {
            for (c = 0; c < numCols; c++) {
               System.out.print (A[l][c] + "  ");
            }
            System.out.println ("");
         }
      System.out.println ("");
   }


   // Left-multiplies lower-triangular matrix Mj by original C_j,
   // where original C_j is in originalMat and result is in genMat.
   // This implementation is safe only if (numRows*(b-1)^2) is a valid int.
   private void leftMultiplyMat (int j, int[][] Mj) {
      int l, c, i, sum;   // Dimension j, row l, column c for new C_j.
      for (l = 0; l < numRows ; l++) {
         for (c = 0; c < numCols; c++) {
            // Multiply row l of M_j by column c of C_j.
            sum = 0;
            for (i = 0; i <= l; i++)
               sum += Mj[l][i] * originalMat[j*numCols+c][i];
            genMat[j*numCols+c][l] = sum % b;
         }
      }
   }


   // Left-multiplies diagonal matrix Mj by original C_j,
   // where original C_j is in originalMat and result is in genMat.
   // This implementation is safe only if (numRows*(b-1)^2) is a valid int.
   private void leftMultiplyMatDiag (int j, int[][] Mj) {
      int l, c, sum;   // Dimension j, row l, column c for new C_j.
      for (l = 0; l < numRows ; l++) {
         for (c = 0; c < numCols; c++) {
            // Multiply row l of M_j by column c of C_j.
            sum = Mj[l][l] * originalMat[j*numCols+c][l];
            genMat[j*numCols+c][l] = sum % b;
         }
      }
   }


   // Right-multiplies original C_j by upper-triangular matrix Mj,
   // where original C_j is in originalMat and result is in genMat.
   // This implementation is safe only if (numCols*(b-1)^2) is a valid int.
   private void rightMultiplyMat (int j, int[][] Mj) {
      int l, c, i, sum;   // Dimension j, row l, column c for new C_j.
      for (l = 0; l < numRows ; l++) {
         for (c = 0; c < numCols; c++) {
            // Multiply row l of C_j by column c of M_j.
            sum = 0;
            for (i = 0; i <= c; i++)
               sum += originalMat[j*numCols+i][l] * Mj[i][c];
            genMat[j*numCols+c][l] = sum % b;
         }
      }
   }


   private int getFaureIndex (String method, int sb, int flag) {
      // Check for errors in ...FaurePermut. Returns the index ib of the
      // base b in primes, i.e.  b = primes[ib].
      if (sb >= b)
         throw new IllegalArgumentException (PrintfFormat.NEWLINE +
            "   sb >= base in " + method);
      if (sb < 1)
         throw new IllegalArgumentException (PrintfFormat.NEWLINE +
            "   sb = 0 in " + method);
      if ((flag > 2) || (flag < 0))
         throw new IllegalArgumentException (
             PrintfFormat.NEWLINE + "   lowerFlag not in {0, 1, 2} in "
             + method);

      // Find index ib of base in array primes
      int ib = 0;
      while ((ib < primes.length) && (primes[ib] < b))
          ib++;
      if (ib >= primes.length)
         throw new IllegalArgumentException ("base too large in " + method);
      if (b != primes[ib])
         throw new IllegalArgumentException (
            "Faure factors are not implemented for this base in " + method);
      return ib;
      }
\end{hide}

   public void leftMatrixScramble (RandomStream stream) \begin{hide} {
      int j, l, c;  // dimension j, row l, column c.

      // If genMat contains the original gen. matrices, copy to originalMat.
      if (originalMat == null) {
         originalMat = genMat;
         genMat = new int[dim * numCols][numRows];
      }

      // Constructs the lower-triangular scrambling matrices M_j, r by r.
      int[][][] scrambleMat = new int[dim][numRows][numRows];
      for (j = 0 ; j < dim; j++) {
         for (l = 0; l < numRows; l++) {
            for (c = 0; c < numRows; c++) {
               if (c == l)                   // No zero on the diagonal.
                  scrambleMat[j][l][c] = stream.nextInt(1, b - 1) ;
               else if (c < l)
                  scrambleMat[j][l][c] = stream.nextInt(0, b - 1);
               else
                  scrambleMat[j][l][c] = 0;  // Zeros above the diagonal;
            }
         }
      }

      // Multiply M_j by the generator matrix C_j for each j.
      for (j = 0 ; j < dim; j++) leftMultiplyMat (j, scrambleMat[j]);
   }\end{hide}
\end{code}
 \begin{tabb}
   Applies a linear scramble by multiplying each $\mathbf{C}_j$ on the left
   by a $w\times w$ nonsingular lower-triangular matrix $\mathbf{M}_j$,
   as suggested by Matou\v{s}ek\latex{ \cite{mMAT99a}} and implemented
   by Hong and Hickernell\latex{ \cite{vHON03a}}. The diagonal entries of
   each matrix $\mathbf{M}_j$ are generated uniformly over
   $\{1,\dots,b-1\}$, the entries below the diagonal are generated uniformly
   over $\{0,\dots,b-1\}$, and all these entries are generated independently.
   This means that in base $b=2$, all diagonal elements are equal to 1.
\richard{Les matrices de \texttt{leftMatrixScramble} sont carr\'ees et
  triangulaires inf\'erieures. PL pense qu'il faut consid\'erer la
  possibilit\'e de rajouter des lignes \`a ces matrices pour pouvoir
  randomiser plus les derniers chiffres ou les derniers bits.}
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used to generate the randomness}
\end{htmlonly}
\begin{code}

   public void leftMatrixScrambleDiag (RandomStream stream) \begin{hide} {
      int j, l, c;  // dimension j, row l, column c.

      // If genMat contains the original gen. matrices, copy to originalMat.
      if (originalMat == null) {
         originalMat = genMat;
         genMat = new int[dim * numCols][numRows];
      }

      // Constructs the diagonal scrambling matrices M_j, r by r.
      int[][][] scrambleMat = new int[dim][numRows][numRows];
      for (j = 0 ; j < dim; j++) {
         for (l = 0; l < numRows; l++) {
            for (c = 0; c < numRows; c++) {
               if (c == l)                   // No zero on the diagonal.
                  scrambleMat[j][l][c] = stream.nextInt(1, b - 1) ;
               else
                  scrambleMat[j][l][c] = 0;  // Diagonal matrix;
            }
         }
      }

      // Multiply M_j by the generator matrix C_j for each j.
      for (j = 0 ; j < dim; j++) leftMultiplyMatDiag (j, scrambleMat[j]);
   }\end{hide}
\end{code}
 \begin{tabb}
   Similar to \method{leftMatrixScramble}{} except that all the
   off-diagonal elements of the $\mathbf{M}_j$ are 0.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used to generate the randomness}
\end{htmlonly}
\begin{code}\begin{hide}

   private void LMSFaurePermut (String method, RandomStream stream, int sb,
      int lowerFlag) {
/*
   If \texttt{lowerFlag = 2}, the off-diagonal elements below the diagonal
   of $\mathbf{M}_j$ are chosen as in \method{leftMatrixScramble}{}.
   If \texttt{lowerFlag = 1}, the off-diagonal elements below the diagonal of
   $\mathbf{M}_j$ are also chosen from the restricted set. If
    \texttt{lowerFlag = 0}, the off-diagonal elements of $\mathbf{M}_j$ are all 0.
*/
      int ib = getFaureIndex (method, sb, lowerFlag);

      // If genMat contains the original gen. matrices, copy to originalMat.
      if (originalMat == null) {
         originalMat = genMat;
         genMat = new int[dim * numCols][numRows];
      }

      // Constructs the lower-triangular scrambling matrices M_j, r by r.
      int jb;
      int j, l, c;  // dimension j, row l, column c.
      int[][][] scrambleMat = new int[dim][numRows][numRows];
      for (j = 0 ; j < dim; j++) {
         for (l = 0; l < numRows; l++) {
            for (c = 0; c < numRows; c++) {
               if (c == l) {
                  jb = stream.nextInt(0, sb - 1);
                  scrambleMat[j][l][c] = FaureFactor[ib][jb];
               } else if (c < l) {
                  if (lowerFlag == 2) {
                     scrambleMat[j][l][c] = stream.nextInt(0, b - 1);
                  } else if (lowerFlag == 1) {
                     jb = stream.nextInt(0, sb - 1);
                     scrambleMat[j][l][c] = FaureFactor[ib][jb];
                  } else {   // lowerFlag == 0
                     scrambleMat[j][l][c] = 0;
                  }
               } else
                  scrambleMat[j][l][c] = 0;  // Zeros above the diagonal;
            }
         }
      }
      // printMat (dim, scrambleMat, method);

      // Multiply M_j by the generator matrix C_j for each j.
      if (lowerFlag == 0)
         for (j = 0 ; j < dim; j++) leftMultiplyMatDiag (j, scrambleMat[j]);
      else
         for (j = 0 ; j < dim; j++) leftMultiplyMat (j, scrambleMat[j]);
   }\end{hide}

   public void leftMatrixScrambleFaurePermut (RandomStream stream, int sb)\begin{hide} {
       LMSFaurePermut ("leftMatrixScrambleFaurePermut", stream, sb, 2);
   }\end{hide}
\end{code}
 \begin{tabb}
   Similar to \method{leftMatrixScramble}{} except that the diagonal elements
   of each matrix $\mathbf{M}_j$ are chosen from a restricted set of the best
   integers as calculated by Faure \cite{rFAU02a}. They are generated
   uniformly over the first \texttt{sb} elements of array $F$, where $F$ is
   made up of a permutation of the integers $[1..(b-1)]$. These integers are
   sorted by increasing  order of the upper bounds of the extreme discrepancy
   for the given integer.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used to generate the randomness}
   \param{sb}{Only the first $sb$ elements of $F$ are used}
\end{htmlonly}
\begin{code}

   public void leftMatrixScrambleFaurePermutDiag (RandomStream stream,
                                                  int sb)\begin{hide} {
       LMSFaurePermut ("leftMatrixScrambleFaurePermutDiag",
                              stream, sb, 0);
   }\end{hide}
\end{code}
 \begin{tabb}
   Similar to \method{leftMatrixScrambleFaurePermut}{} except that all
   off-diagonal elements are 0.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used to generate the randomness}
   \param{sb}{Only the first $sb$ elements of $F$ are used}
\end{htmlonly}
\begin{code}

   public void leftMatrixScrambleFaurePermutAll (RandomStream stream,
                                                 int sb)\begin{hide} {
       LMSFaurePermut ("leftMatrixScrambleFaurePermutAll",
                              stream, sb, 1);
   }\end{hide}
\end{code}
 \begin{tabb}
   Similar to \method{leftMatrixScrambleFaurePermut}{} except that the
   elements under the diagonal are also
   chosen from the same restricted set as the diagonal elements.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used to generate the randomness}
   \param{sb}{Only the first $sb$ elements of $F$ are used}
\end{htmlonly}
\begin{code}

   public void iBinomialMatrixScramble (RandomStream stream)\begin{hide} {
      int j, l, c;  // dimension j, row l, column c.
      int diag;     // random entry on the diagonal;
      int col1;     // random entries below the diagonal;

      // If genMat is original generator matrices, copy it to originalMat.
      if (originalMat == null) {
         originalMat = genMat;
         genMat = new int[dim * numCols][numRows];
      }

      // Constructs the lower-triangular scrambling matrices M_j, r by r.
      int[][][] scrambleMat = new int[dim][numRows][numRows];
      for (j = 0 ; j < dim; j++) {
         diag = stream.nextInt(1, b - 1);
         for (l = 0; l < numRows; l++) {
            // Single random nonzero element on the diagonal.
            scrambleMat[j][l][l] = diag;
            for (c = l+1; c < numRows; c++) scrambleMat[j][l][c] = 0;
         }
         for (l = 1; l < numRows; l++) {
            col1 = stream.nextInt(0, b - 1);
            for (c = 0; l+c < numRows; c++) scrambleMat[j][l+c][c] = col1;
         }
      }
      // printMat (dim, scrambleMat,  "iBinomialMatrixScramble");
      for (j = 0 ; j < dim; j++) leftMultiplyMat (j, scrambleMat[j]);
   }\end{hide}
\end{code}
 \begin{tabb}
   Applies the $i$-binomial matrix scramble proposed by Tezuka \cite{rTEZ02a}
   \latex{ (see also \cite{vOWE03a})}.
   This multiplies each $\mathbf{C}_j$ on the left
   by a $w\times w$ nonsingular lower-triangular matrix $\mathbf{M}_j$ as in
   \method{leftMatrixScramble}{}, but with the additional constraint that
   all entries on any given diagonal or subdiagonal of $\mathbf{M}_j$ are identical.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used as generator of the randomness}
\end{htmlonly}
\begin{code}\begin{hide}

   private void iBMSFaurePermut (String method, RandomStream stream,
                                 int sb, int lowerFlag) {
      int ib = getFaureIndex (method, sb, lowerFlag);

      // If genMat is original generator matrices, copy it to originalMat.
      if (originalMat == null) {
         originalMat = genMat;
         genMat = new int[dim * numCols][numRows];
      }

      // Constructs the lower-triangular scrambling matrices M_j, r by r.
      int j, l, c;  // dimension j, row l, column c.
      int diag;     // random entry on the diagonal;
      int col1;     // random entries below the diagonal;
      int jb;
      int[][][] scrambleMat = new int[dim][numRows][numRows];
      for (j = 0 ; j < dim; j++) {
         jb = stream.nextInt(0, sb - 1);
         diag = FaureFactor[ib][jb];
         for (l = 0; l < numRows; l++) {
            // Single random nonzero element on the diagonal.
            scrambleMat[j][l][l] = diag;
            for (c = l+1; c < numRows; c++) scrambleMat[j][l][c] = 0;
         }
         for (l = 1; l < numRows; l++) {
            if (lowerFlag == 2) {
               col1 = stream.nextInt(0, b - 1);
            } else if (lowerFlag == 1) {
               jb = stream.nextInt(0, sb - 1);
               col1 = FaureFactor[ib][jb];
            } else {   // lowerFlag == 0
               col1 = 0;
            }
            for (c = 0; l+c < numRows; c++) scrambleMat[j][l+c][c] = col1;
         }
      }
      // printMat (dim, scrambleMat, method);

      if (lowerFlag > 0)
         for (j = 0 ; j < dim; j++) leftMultiplyMat (j, scrambleMat[j]);
      else
         for (j = 0 ; j < dim; j++) leftMultiplyMatDiag (j, scrambleMat[j]);
   }\end{hide}

   public void iBinomialMatrixScrambleFaurePermut (RandomStream stream,
                                                   int sb)\begin{hide} {
       iBMSFaurePermut ("iBinomialMatrixScrambleFaurePermut",
                              stream, sb, 2);
   }\end{hide}
\end{code}
 \begin{tabb}
   Similar to \method{iBinomialMatrixScramble}{} except that the diagonal
   elements of each matrix $\mathbf{M}_j$ are chosen as in
  \method{leftMatrixScrambleFaurePermut}{}.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used to generate the randomness}
   \param{sb}{Only the first $sb$ elements of $F$ are used}
\end{htmlonly}
\begin{code}

   public void iBinomialMatrixScrambleFaurePermutDiag (RandomStream stream,
                                                       int sb)\begin{hide} {
       iBMSFaurePermut ("iBinomialMatrixScrambleFaurePermutDiag",
                               stream, sb, 0);
   }\end{hide}
\end{code}
 \begin{tabb}
   Similar to \method{iBinomialMatrixScrambleFaurePermut}{} except that all the
   off-diagonal elements are 0.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used to generate the randomness}
   \param{sb}{Only the first $sb$ elements of $F$ are used}
\end{htmlonly}
\begin{code}

   public void iBinomialMatrixScrambleFaurePermutAll (RandomStream stream,
                                                      int sb)\begin{hide} {
       iBMSFaurePermut ("iBinomialMatrixScrambleFaurePermutAll",
                               stream, sb, 1);
   }\end{hide}
\end{code}
 \begin{tabb}
   Similar to \method{iBinomialMatrixScrambleFaurePermut}{} except that the
   elements under the diagonal are also
   chosen from the same restricted set as the diagonal elements.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used to generate the randomness}
   \param{sb}{Only the first $sb$ elements of $F$ are used}
\end{htmlonly}
\begin{code}

   public void stripedMatrixScramble (RandomStream stream)\begin{hide} {
      int j, l, c;  // dimension j, row l, column c.
      int diag;     // random entry on the diagonal;
      int col1;     // random entries below the diagonal;

      // If genMat contains original gener. matrices, copy it to originalMat.
      if (originalMat == null) {
         originalMat = genMat;
         genMat = new int[dim * numCols][numRows];
      }

      // Constructs the lower-triangular scrambling matrices M_j, r by r.
      int[][][] scrambleMat = new int[dim][numRows][numRows];
      for (j = 0 ; j < dim; j++) {
         for (c = 0; c < numRows; c++) {
            diag = stream.nextInt (1, b - 1);   // Random entry in this column.
            for (l = 0; l < c; l++)         scrambleMat[j][l][c] = 0;
            for (l = c; l < numRows; l++) scrambleMat[j][l][c] = diag;
         }
      }
      // printMat (dim, scrambleMat,  "stripedMatrixScramble");
      for (j = 0 ; j < dim; j++) leftMultiplyMat (j, scrambleMat[j]);
   }\end{hide}
\end{code}
 \begin{tabb}
   Applies the \emph{striped matrix scramble} proposed by Owen \cite{vOWE03a}.
   It multiplies each $\mathbf{C}_j$ on the left
   by a $w\times w$ nonsingular lower-triangular matrix $\mathbf{M}_j$ as in
   \method{leftMatrixScramble}{}, but with the additional constraint that
   in any column, all entries below the diagonal are equal to the
   diagonal entry, which is generated randomly over $\{1,\dots,b-1\}$.
   Note that for $b=2$, the matrices $\mathbf{M}_j$ become deterministic, with
   all entries on and below the diagonal equal to 1.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used as generator of the randomness}
\end{htmlonly}
\begin{code}

   public void stripedMatrixScrambleFaurePermutAll (RandomStream stream,
                                                    int sb)\begin{hide} {
      int ib = getFaureIndex ("stripedMatrixScrambleFaurePermutAll", sb, 1);

      // If genMat contains original gener. matrices, copy it to originalMat.
      if (originalMat == null) {
         originalMat = genMat;
         genMat = new int[dim * numCols][numRows];
      }

      // Constructs the lower-triangular scrambling matrices M_j, r by r.
      int j, l, c;  // dimension j, row l, column c.
      int diag;     // random entry on the diagonal;
      int col1;     // random entries below the diagonal;
      int jb;
      int[][][] scrambleMat = new int[dim][numRows][numRows];
      for (j = 0 ; j < dim; j++) {
         for (c = 0; c < numRows; c++) {
            jb = stream.nextInt(0, sb - 1);
            diag = FaureFactor[ib][jb];  // Random entry in this column.
            for (l = 0; l < c; l++)         scrambleMat[j][l][c] = 0;
            for (l = c; l < numRows; l++) scrambleMat[j][l][c] = diag;
         }
      }
      // printMat (dim, scrambleMat,  "stripedMatrixScrambleFaurePermutAll");
      for (j = 0 ; j < dim; j++) leftMultiplyMat (j, scrambleMat[j]);
   }\end{hide}
\end{code}
 \begin{tabb}
   Similar to \method{stripedMatrixScramble}{} except that the
   elements on and under the diagonal of each matrix $\mathbf{M}_j$ are
   chosen as in \method{leftMatrixScrambleFaurePermut}{}.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used as generator of the randomness}
   \param{sb}{Only the first $sb$ elements of $F$ are used}
\end{htmlonly}
\begin{code}

   public void rightMatrixScramble (RandomStream stream) \begin{hide} {
      int j, c, l, i, sum;  // dimension j, row l, column c, of genMat.

      // SaveOriginalMat();
      if (originalMat == null) {
         originalMat = genMat;
         genMat = new int[dim * numCols][numRows];
      }

      // Generate an upper triangular matrix for Faure-Tezuka right-scramble.
      // Entry [l][c] is in row l, col c.
      int[][] scrambleMat = new int[numCols][numCols];
      for (c = 0; c < numCols; c++) {
         for (l = 0; l < c; l++) scrambleMat[l][c] = stream.nextInt (0, b - 1);
         scrambleMat[c][c] = stream.nextInt (1, b - 1);
         for (l = c+1; l < numCols; l++) scrambleMat[l][c] = 0;
      }

      // printMat0 (scrambleMat,  "rightMatrixScramble");
      // Right-multiply the generator matrices by the scrambling matrix.
      for (j = 0 ; j < dim; j++) rightMultiplyMat (j, scrambleMat);
   }\end{hide}
\end{code}
 \begin{tabb}
   Applies a linear scramble by multiplying each $\mathbf{C}_j$ on the right
   by a single $k\times k$ nonsingular upper-triangular matrix $\mathbf{M}$,
   as suggested by Faure and Tezuka\latex{ \cite{rFAU02a} (see also
   \cite{vHON03a})}.
   The diagonal entries of the matrix $\mathbf{M}$ are generated uniformly over
   $\{1,\dots,b-1\}$, the entries above the diagonal are generated uniformly
   over $\{0,\dots,b-1\}$, and all the entries are generated independently.
   The effect of this scramble is only to change the order in which the
   points are generated.  If one computes the average value of a
   function over \emph{all} the points of a given digital net, or over
   a number of points that is a power of the basis, then this
   scramble makes no difference.
\end{tabb}
\begin{htmlonly}
   \param{stream}{random number stream used as generator of the randomness}
\end{htmlonly}
\begin{hide} %%%%
\begin{code}

   public void unrandomize() {
      resetGeneratorMatrices();
      digitalShift = null;
   }
\end{code}
\begin{tabb}
   Restores the original generator matrices and
   removes the random shift.
\end{tabb}
\end{hide}
\begin{code}

   public void resetGeneratorMatrices()\begin{hide} {
      if (originalMat != null) {
         genMat = originalMat;
         originalMat = null;
      }
   }\end{hide}
\end{code}
\begin{tabb}
   Restores the original generator matrices.
   This removes the current linear matrix scrambles.
\end{tabb}
\begin{code}

   public void eraseOriginalGeneratorMatrices()\begin{hide} {
      originalMat = null;
   }\end{hide}
\end{code}
\begin{tabb}
   Erases the original generator matrices and replaces them by
   the current ones.  The current linear matrix scrambles thus become
   \emph{permanent}.  This is useful if we want to apply several
   scrambles in succession to a given digital net.
\end{tabb}
\begin{code}

   public void printGeneratorMatrices (int s) \begin{hide} {
      // row l, column c, dimension j.
      for (int j = 0; j < s; j++) {
         System.out.println ("dim = " + (j+1) + PrintfFormat.NEWLINE);
         for (int l = 0; l < numRows; l++) {
            for (int c = 0; c < numCols; c++)
               System.out.print (genMat[j*numCols+c][l] + "  ");
            System.out.println ("");
         }
         System.out.println ("----------------------------------");
      }
   }\end{hide}
\end{code}
\begin{tabb}
  Prints the generator matrices in standard form for dimensions 1 to $s$.
\end{tabb}

\begin{code}\begin{hide}

      // Computes the digital expansion of $i$ in base $b$, and the digits
      // of the Gray code of $i$.
      // These digits are placed in arrays \texttt{bary} and \texttt{gray},
      // respectively, and the method returns the number of digits in the
      // expansion.  The two arrays are assumed to be of sizes larger or
      // equal to this new number of digits.
      protected int intToDigitsGray (int b, int i, int numDigits,
                                            int[] bary, int[] gray) {
         if (i == 0)
            return 0;
         int idigits = 0; // Num. of digits in b-ary and Gray repres.
         int c;
         // convert i to b-ary representation, put digits in bary[].
         for (c = 0; i > 0; c++) {
            idigits++;
            bary[c] = i % b;
            i = i / b;
         }
         // convert b-ary representation to Gray code.
         gray[idigits-1] = bary[idigits-1];
         int diff;
         for (c = 0;  c < idigits-1; c++) {
            diff = bary[c] - bary[c+1];
            if (diff < 0)
               gray[c] = diff + b;
            else
               gray[c] = diff;
         }
         for (c = idigits; c < numDigits; c++) gray[c] = bary[c] = 0;
         return idigits;
      }


// ************************************************************************

   protected class DigitalNetIterator extends
          PointSet.DefaultPointSetIterator {

      protected int idigits;        // Num. of digits in current code.
      protected int[] bdigit;       // b-ary code of current point index.
      protected int[] gdigit;       // Gray code of current point index.
      protected int dimS;           // = dim except in the shift iterator
                                    // children where it is = dim + 1.
      protected int[] cachedCurPoint; // Digits of coords of the current point,
                           // with digital shift already applied.
                           // Digit l of coord j is at pos. [j*outDigits + l].
                           // Has one more dimension for the shift iterators.

      public DigitalNetIterator() {
 //        EpsilonHalf = 0.5 / Math.pow ((double) b, (double) outDigits);
         bdigit = new int[numCols];
         gdigit = new int[numCols];
         dimS = dim;
         cachedCurPoint = new int[(dim + 1)* outDigits];
         init();  // This call is important so that subclasses don't
                  // automatically call 'resetCurPointIndex' at the time
                  // of construction as this may cause a subtle
                  // 'NullPointerException'
      }

      public void init() { // See constructor
         resetCurPointIndex();
      }

      // We want to avoid generating 0 or 1
      public double nextDouble() {
         return nextCoordinate() + EpsilonHalf;
      }

      public double nextCoordinate() {
         if (curPointIndex >= numPoints || curCoordIndex >= dimS)
            outOfBounds();
         int start = outDigits * curCoordIndex++;
         double sum = 0;
         // Can always have up to outDigits digits, because of digital shift.
         for (int k = 0; k < outDigits; k++)
            sum += cachedCurPoint [start+k] * factor[k];
         if (digitalShift != null)
            sum += EpsilonHalf;
        return sum;
      }

      public void resetCurPointIndex() {
         if (digitalShift == null)
            for (int i = 0; i < cachedCurPoint.length; i++)
               cachedCurPoint[i] = 0;
         else {
            if (dimShift < dimS)
               addRandomShift (dimShift, dimS, shiftStream);
            for (int j = 0; j < dimS; j++)
               for (int k = 0; k < outDigits; k++)
                  cachedCurPoint[j*outDigits + k] = digitalShift[j][k];
         }
         for (int i = 0; i < numCols; i++) bdigit[i] = 0;
         for (int i = 0; i < numCols; i++) gdigit[i] = 0;
         curPointIndex = 0;
         curCoordIndex = 0;
         idigits = 0;
      }

      public void setCurPointIndex (int i) {
         if (i == 0) {
            resetCurPointIndex();
            return;
         }
         curPointIndex = i;
         curCoordIndex = 0;
         if (digitalShift != null && dimShift < dimS)
             addRandomShift (dimShift, dimS, shiftStream);

         // Digits of Gray code, used to reconstruct cachedCurPoint.
         idigits = intToDigitsGray (b, i, numCols, bdigit, gdigit);
         int c, j, l, sum;
         for (j = 0; j < dimS; j++) {
            for (l = 0; l < outDigits; l++) {
               if (digitalShift == null)
                  sum = 0;
               else
                  sum = digitalShift[j][l];
               if (l < numRows)
                  for (c = 0; c < idigits; c++)
                     sum += genMat[j*numCols+c][l] * gdigit[c];
               cachedCurPoint [j*outDigits+l] = sum % b;
            }
         }
      }

      public int resetToNextPoint() {
         // incremental computation.
         curPointIndex++;
         curCoordIndex = 0;
         if (curPointIndex >= numPoints)
            return curPointIndex;

         // Update the digital expansion of i in base b, and find the
         // position of change in the Gray code. Set all digits == b-1 to 0
         // and increase the first one after by 1.
         int pos;      // Position of change in the Gray code.
         for (pos = 0; gdigit[pos] == b-1; pos++)
            gdigit[pos] = 0;
         gdigit[pos]++;

         // Update the cachedCurPoint by adding the column of the gener.
         // matrix that corresponds to the Gray code digit that has changed.
         // The digital shift is already incorporated in the cached point.
         int c, j, l;
         int lsup = numRows;        // Max index l
         if (outDigits < numRows)
            lsup = outDigits;
         for (j = 0; j < dimS; j++) {
            for (l = 0; l < lsup; l++) {
               cachedCurPoint[j*outDigits + l] += genMat[j*numCols + pos][l];
               cachedCurPoint[j * outDigits + l] %= b;
            }
         }
         return curPointIndex;
      }
   }


// ************************************************************************

   protected class DigitalNetIteratorNoGray extends DigitalNetIterator {

      public DigitalNetIteratorNoGray() {
         super();
      }

      public void setCurPointIndex (int i) {
         if (i == 0) {
            resetCurPointIndex();
            return;
         }
         curPointIndex = i;
         curCoordIndex = 0;
         if (dimShift < dimS)
             addRandomShift (dimShift, dimS, shiftStream);

         // Convert i to b-ary representation, put digits in bdigit.
         idigits = intToDigitsGray (b, i, numCols, bdigit, gdigit);
         int c, j, l, sum;
         for (j = 0; j < dimS; j++) {
            for (l = 0; l < outDigits; l++) {
               if (digitalShift == null)
                  sum = 0;
               else
                  sum = digitalShift[j][l];
               if (l < numRows)
                  for (c = 0; c < idigits; c++)
                     sum += genMat[j*numCols+c][l] * bdigit[c];
               cachedCurPoint [j*outDigits+l] = sum % b;
            }
         }
      }

      public int resetToNextPoint() {
         curPointIndex++;
         curCoordIndex = 0;
         if (curPointIndex >= numPoints)
            return curPointIndex;

         // Find the position of change in the digits of curPointIndex in base
         // b. Set all digits = b-1 to 0; increase the first one after by 1.
         int pos;
         for (pos = 0; bdigit[pos] == b-1; pos++)
            bdigit[pos] = 0;
         bdigit[pos]++;

         // Update the digital expansion of curPointIndex in base b.
         // Update the cachedCurPoint by adding 1 unit at the digit pos.
         // If pos > 0, remove b-1 units in the positions < pos. Since
         // calculations are mod b, this is equivalent to adding 1 unit.
         // The digital shift is already incorporated in the cached point.
         int c, j, l;
         int lsup = numRows;        // Max index l
         if (outDigits < numRows)
            lsup = outDigits;
         for (j = 0; j < dimS; j++) {
            for (l = 0; l < lsup; l++) {
               for (c = 0; c <= pos; c++)
                  cachedCurPoint[j*outDigits + l] += genMat[j*numCols + c][l];
               cachedCurPoint[j * outDigits + l] %= b;
            }
         }
         return curPointIndex;
      }
   }

}\end{hide}
\end{code}
